%Richtmyer Scheme
function u=Richt(u0,dt,maxa,N,al,b)
u=zeros(1,N);
h=1/N;
la=dt/h;
flux = @(uu) uu.^al/b;
for j=1:N 
    if j~=1&&j~=N
        ustar=(u0(j+1)+u0(j))/2-la*(flux(u0(j+1))-flux(u0(j)))/2;
        fp=flux(ustar);
        
        ustar=(u0(j)+u0(j-1))/2-la*(flux(u0(j))-flux(u0(j-1)))/2;
        fn=flux(ustar);
    elseif j==1
        ustar=(u0(j+1)+u0(j))/2-la*(flux(u0(j+1))-flux(u0(j)))/2;
        fp=flux(ustar);
        
        ustar=(u0(j)+u0(N))/2-la*(flux(u0(j))-flux(u0(N)))/2;
        fn=flux(ustar);
    else
        ustar=(u0(1)+u0(j))/2-la*(flux(u0(1))-flux(u0(j)))/2;
        fp=flux(ustar);
        
        ustar=(u0(j)+u0(j-1))/2-la*(flux(u0(j))-flux(u0(j-1)))/2;
        fn=flux(ustar);
    end
    u(j)=u0(j)-la*(fp-fn);
end